Solving for pressure : acceleration and accuracy

Meso and large scale models tend nowadays to be non-hydrostatic. Time step restrictions due to fast waves are alleviated by considering semi-implicit schemes (Thomas et al, 2003) or by use of the anelastic approximation (Meso-NH model, Lafore et al, 1998, or Clark, 1977 among others); in both cases an elliptic equation has to be solved. This equation is typically ill conditioned, due to the range of scales considered and the ratio of horizontal to vertical scales; it is also non-separable, as the domain is not parallepipedic, due to the presence of orography; in most discretizations it is also found non-symmetric. Thomas et al (2003) make a survey of the different methods to solve this equation; they adopt the preconditioning proposed in Bernardet (1995), hereafter B95 and rely on generalizations of the conjugate gradient method to non-symmetric problems such as GCR(k). In B95 it was noticed that the divergence in the continuity equation and the gradient of pressure are naturally adjoint operators, but, if the divergence is naturally discretized in flux form, two discretizations of the pressure gradient come naturally, and moreover the elliptic equation might be discretized directly without reference to the gradient or divergence operators, with a minimum use of averaging for example. With a certain placement of the quantities defining the metrics of the computational grid, it was found that the discretizations of the pressure terms are equivalent (outside boundaries), but that disposition was not used in the Clark (1977) model nor retained in the Meso-NH model. We therefore seek in this paper to convince the reader that a symmetric elliptic equation can be obtained by an adequate design of the extrapolations at the boundary for the discretized gradient and divergence operators; the Helmholtz equation will then be solved at a minimal cost; for example, the orthomin solver is adequate for non-symmetric problems but necessitates two applications of the preconditioner per iteration, instead of one for the standard conjugate gradient.